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ABSTRACT 


Some  aspects  of  the  numerical  solution  of  the  Tricomi 

equat i on 

y  $  -  $  =*0 

xx  yy 

and  the  inverted  Tricomi  equation 

y$  -  $  *  0 

yy  xx 

with  particular  emphasis  on  periodic  problems  are  studied. 


(  1  ) 


(2) 
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SIGNIFICANCE  AND  EXPLANATION 
We  consider  in  this  paper  some  aspects  of  the  numerical 
solution  of  the  Tricomi  equation 

y$  -  <P  =  0  (  1  ) 

xx  yy 

and  the  inverted  Tricomi  equation 

y<J>  -  <J>  =  o  (2) 

yy  xx 

with  particular  emphasis  on  periodic  problems.  These  periodic 
problems  are  of  actual  physical  interest:  the  former  is  a  model 
for  the  deflection  of  a  floppy  disc  considered  as  a  rotating 
membrane,  while  the  latter  is  a  model  for  the  transonic  deLaval 
nozzle.  Since  most  studies  of  the  Tricomi  equations  have  been  in 
domains  bounded  by  one  or  more  characteristics,  such  periodic 
problems  offer  some  different  viewpoints  and  some  different 
qualitative  insight  into  these  mixed  elliptic-hyperbolic 
equat i ons . 

We  shall  test  out  various  numerical  schemes  on  these  two 
problems.  These  equations  being  linear,  many  mathematical 
questions,  e.g.  unique  solvability,  convergence,  etc.,  can  be 
easily  answered  on  these  model  problems  that  would  help  shed 
light  on  actual  nonlinear  numerical  procedures,  say,  those  used 
L..nsonic  flow.  We  shall  concentrate  on  the  unique 
solvability  of  the  algebraic  equations,  and  expect  convergence 
,-tc.  can  be  studied  equally  readily. 


The  responsibility  for  the  wording  and  views  expressed  in  this 
descriptive  summary  lies  with  MRC ,  and  not  with  the  authors  of 
this  report. 


SOME  REMARKS  ON  THE  NUMERICAL  SOLUTION  OP  TRICOMI-TYPE  EQUATIONS 
C.  K.  Chu*,  L.  W.  Xiang  and  Z.  K.  Yao 

1.  INTRODUCTION 

We  consider  in  this  paper  some  aspects  of  the  numerical 
solution  of  the  Tricomi  equation 

yd>  -  <j>  =0  (1) 

,yxx  Yyy 

and  the  inverted  Tricomi  equation 

yd)  -  $  =  0  (2) 

with  particular  emphasis  on  periodic  problems.  These  peri¬ 
odic  problems  are  of  actual  physical  interest :  the  former 
is  a  model  for  the  deflection  of  a  floppy  disc  considered  as 
a  rotating  membrane,  while  the  latter  is  a  model  for  the 
transonic  deLaval  nozzle.  Since  most  studies  of  the  Tricomi 
equation  have  been  in  domains  bounded  by  one  or  more  charac¬ 
teristics,  such  periodic  problems  offer  some  different  view¬ 
points  and  some  different  qualitative  insight  into  these 
mixed  elliptic-hyperbolic  equations. 

We  shall  test  out  various  numerical  schemes  on  these 
two  problems.  These  equations  being  linear,  many  mathema¬ 
tical  questions,  e.g.  unique  solvability,  convergence,  etc., 
can  be  easily  answered  on  these  model  problems  that  would 
help  shed  light  on  actual  nonlinear  numerical  procedures, 
say,  those  used  in  transonic  flow.  We  shall  concentrate  on 
the  unique  solvability  of  the  algebraic  equations,  and  expect 
convergence  etc.  can  be  studied  equally  readily. 

fPartially  supported  by  USDOE  Contract  No.  DE-AC  02-76ET53016  at 
Columbia  University. 

Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80- 
C-0041  . 


2.  THE  PERIODIC  TRICOMI  PROBLEM  AND  FLOPPY  DISC 

A  rotating  circular  membrane  clamped  at  the  inner  edge 
r  *  a  and  loaded  with  a  transverse  force  (fig.  1)  is  des¬ 
cribed  by  the  following  equation  |l| 

-  (1-r2)  wrr  -  i  (3r2-  1)  wr  +  \2  (3r2-  1)  wQe  =  F  (3) 

Here  w  is  the  (small)  deflection  transverse  to  the  disc, 

0  is  the  usual  angular  coordinate,  r  the  radial  coordinate 
normalized  with  respect  to  the  outer  radius,  and  the  driving 
term  F  contains  the  loading  force  per  unit  area,  angular 
velocity,  disc  thickness,  and  material  constants.  The 
appropriate  boundary  conditions  are 

w  (a, 8)  =  0 

w  (1, 6)  <  *  (4) 

w(r,8)  =  w(r,  0  +  2  tj) 

expressing,  respectively,  clamping  at  the  inner  edge,  finite 
deflection  at  the  outer  edge,  and  periodicity  in  0.  The 
equation  is  very  interesting  in  that  it  is  elliptic  in  the 
region  a  <  r  <l//3  ,  hyperbolic  in  the  region  1//3  <  r  <  1, 
and  degenerate  parabolic  on  the  transition  line  r  =  1//3 
and  again  on  the  outer  edge  r  =  1.  The  characteristics  in 
the  hyperbolic  region  have  the  generai  shapes  shown  in  fig.  1. 


Fig.  1.  (Left)  Geometry  of  floppy  disc  problem 
(Right)  Model  using  Tricomi  equation 
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This  problem  can  obviously  be  modeled  by  the  following 
periodic  problem  for  the  Tricomi  equation:  4>  is  to  satisfy 
equation  (1)  in  the  domain  -*>  <x  <  +°°  ,  -a  <  y  <  +b  and  the 
boundary  conditions 

<j>  (or  <p  )  prescribed  at  y  =  -a,  y  =  +b 

y  (5) 

4>  (x,y)  =  4>(x+2tt  ,  y) 


The  modeling  is  not  exact  in  that  the  inessential  first  deri¬ 
vative  term  in  (3)  is  discarded  and  that  the  outer  boundary 
y  =  +b  is  not  a  degenerate  parabolic  line  any  more.  Other 
than  that,  the  essential  features,  including  the  shape  of 
the  characteristics,  are  all  properly  modeled.  The  inhomo¬ 
geneous  differential  equation  (3)  with  homogeneous  boundary 
conditions  (4)  has  also  been  changed  to  the  homogeneous 
equation  (1)  and  inhomogeneous  boundary  conditions  (5),  the 
equivalence  of  these  problems  being  well  known. 

The  problem  is  obviously  treated  by  separation  of  vari¬ 
ables.  Letting 

.  _  ,  inx 

<P  <j>n  e 

the  Fourier  component  4>n  for  each  wave  number  n  (n>  0) 
satisfies  an  Airy  equation 

<J>n"  +  n2y  <J>n  =  0  (6> 

The  solution  to  (6) is  given  in  terms  of  the  Airy  functions 
4>n  =  Aq  Ai(z)  +  Bn  Bi( z)  ,  (7) 

where  z  =  n_2^3y,  and  Afl ,  Bn  are  arbitrary  constants.  For 
n  *  0,  4>  *  A^  +  B_v  trivially.  For  large  z,  the  functions 


^  Ai( z)  and  Bi(z)  behave  like  slightly  damped  exponentials 
or  slightly  damped  sinusoids,  depending  on  whether  z  is 
positive  or  negative: 

AiU)  *  &  e-W3>*3/2  .  *  ..  « 


sin  li  (-z)3/3 


z  <  0 


Bi(«)  -A  *'l/4  eH2'3)* 

■v  cos  (§  <-z>3/3 


z  >  0 


-3- 


The  boundary  conditions  to  be  satisfied  by  ^  are  the  n-th 
Fourier  coefficients  of  the  given  data  at  y  =  -a  and  +b. 
Thus,  the  problem  is  uniquely  solvable  if  and  only  if  the 
determinant 

,2/3 


Ai(n 


For  large 


Ai(n2/3a) 


-2/3 


b) 


Bi( n  '  a) 
Bi(n-2/3b) 


(9) 


n 


1  (4/3)na 

2  e 


this  is  approximately  equivalent  to 
3/2 

COt  ( g  HU  ^ 


I  nb3/2  +  l  )  f  1 


For  given  domains,  condition  (9)  will  be  violated  for 
some  very  large  n.  Thus  one  may  conclude  that  the  floppy 
disc  problem  admits  no  unique  solution.  However,  the  mem¬ 
brane  approximation  in  equation  (3)  is  really  only  valid  for 
n  not  too  large,  say,  n  <  N;  since  if  the  azimuthal  wave 
lengths  become  too  small  compared  to  the  radial  dimension 
of  the  disc,  the  membrane  approximation  fails  as  bending  must 
now  be  taken  into  account.  With  the  limitation  of  n  ^  N  , 
then  the  periodic  problem  in  general  has  a  unique  solution 
provided  (9)  is  satisfied.  For  example,  if  a  »  b  *  1  , 
actual  calculation  shows  that  (9)  is  still  satisfied  for  n 
as  high  as  100. 


3.  THE  INVERTED  TRICOMI  EQUATION  AND  THE  NOZZLE  PROBLEM 

The  converging-diverging  nozzle,  in  which  a  gas  expands 
from  subsonic  velocity  to  supersonic  velocity,  can  obviously 
be  viewed  also  as  a  periodic  problem,  if  we  introduce  poten¬ 
tial  or  streamline  coordinates,  in  which  the  walls  would  be 
either  constant  function  values  or  zero  normal  derivatives 
and  thus  the  function  can  be  extended  to  a  periodic  function 
in  the  usual  elementary  way.  If,  as  is  often  tempting  to  do, 
we  try  to  think  of  modeling  the  nozzle  flow  by  the  Tricomi 
equation,  then  the  arguments  of  the  previous  section  would 
convince  us  that  two  boundary  conditions  are  needed,  one  at 
the  subsonic  upstream  and  one  at  the  supersonic  downstream. 
But  in  fact,  this  is  never  the  case  in  practice:  a  supersonic 
boundary  condition  is  never  prescribed  (see  e.g.  |2|).  The 
'physical"argments  that  signals  do  not  propagate  upstream  are 
inadequate  to  explain  this  discrepancy. 
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The  real  reason  is  that  the  Trieomi  equation  is  the 
wrong  model  to  use  i'or  transonic  flow  in  the  physical  plane, 
although  it  is  the  correct  model  in  the  hodograph  plane.  The 
characteristics  in  the  supersonic  flow  region,  fig.  2,  arc 
perpendicular  to  the  streamlines  at  the  sonic  transition,  and 
therefore,  they  are  almost  tangent  to  the  transition  line. 

The  characteristics  of  the  Trieomi  equation,  on  the  contrary, 
are  perpendicular  to  the  transition  line,  as  in  fig.  1.  Thus 
for  calculating  in  the  physical  plane,  the  Trieomi  equation 
does  not  model  the  nozzle  flow  problem  correctly.  In  fact, 
the  inverted  Trieomi  equation  (2)  is  the  right  model  to  study. 
In  his  book  |3|  ,  Trieomi  also  considered  the  inverted  Tri- 
comi  equation  but  stated  that  the  equation  did  not  seem  to 
have  a  physical  application.  The  present  model  appears  to 
be  indeed  such  an  application. 

The  proper  periodic  boundary  value  problem  for  the 
inverted  Trieomi  equation  (2)  is:  $  is  to  satisfy  (2)  in 
the  domain  -®<  x  <  «°,  -a  <  y  <  b  and  the  boundary  condi¬ 


tions 


«j>  or  <)>y  prescribed  at  y  ■  -a 
No  data  to  be  prescribed  at  y  *  +b 

♦  (x.y)  =  <|>(x+2n .  y) 


(10) 


Again  letting  $  *  <j>ne  ,  the  equation  satisfied  by  each 
is  now 


+  n2 


(11) 


In  the  neighborhood  of  v  *  0,  elementary  series  methods  give 

two  linearly  independent  solutions  for  each  n  f  0: 

„  (1)  «  v  M  n£  v  +  n!  2  _n£_  3  * 

*n  y  21  y  +  312 .’  y  4131  y  +  2) 

/  O  v  /  1  \  V  / 


log  y 


Since  the  derivatives  of  model  velocity  components,  we  re- 

(2) 

quire  that  they  be  square- integrable .  The  solution  <f>n v 
fails  to  satisfy  this  requirement,  and  must  be  rejected.  We 
are  then  left  with  one  solution  in  (12),  corresponding  to 
the  one  boundary  condition  at  y  *  -a  and  no  boundary  condi¬ 
tion  at  y  *  +b.  That  the  converse  cannot  be  prescribed  is 


-5- 


based  on  the  ellipticity  or  hyperbolicity  of  the  original 
partial  differential  equation,  and  cannot  be  seen  from  such 
a  simple  argument  as  just  presented. 

On  the  other  hand,  for  n  =  0,  we  have  <j>Q  ®  +  c^y. 

C 2  is  determined  from  the  Fourier  coeffients  of  the  given 
data  at  y  =  -a,  just  as  all  the  coefficients  of  the  $  ,  n^O, 
are  determined.  is  determined  by  prescribing  the  average 

value  of  <}>(x,0)  on  the  transition  line,  where  <j>  =  o  and 

<J>  is  a  linear  function  in  x. 


Fig.  2.  Sketch  of  deLaval  nozzle  showing 

transition  line  and  characteristics. 


Fig.  3.  Domain  and  grid  for  inverted 
Tricomi  equation  in  model  nozzle  prob¬ 
lem. 
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4.  FINITE  DIFFERENCE  SOLUTION  OF  THE  MODEL  NOZZLE  PROBLEM 


We  apply  the  method  of  Murman  and  Cole  |4|  to  the  model 
nozzle  problem  with  the  inverted  Tricomi  equation.  We  lay 
down  a  finite  difference  grid  (fig.  3)  of  grid  width  Ax  and 
Ay  respectively,  in  which  the  line  y  =  0  is  a  row  of  grid 
points,  to  start;  later  on  this  inessential  restriction  will 
be  removed.  The  Murman-Cole  method  uses  five-point  centered 
differences  if  the  center  point  lies  in  the  elliptic  region 
(as  in  a  Laplace  equation),  and  uses  backward  second  differ¬ 
ences  in  y  if  the  point  in  question  lies  in  the  hyperbolic 
region  (as  in  an  implicit  scheme  lor  the  wave  equation). 

Thus,  we  have 


We  shall  prescribe  boundary  conditions  at  y  -  -a,  and 
prescribe  periodic  boundary  conditions  on  the  side  walls. 

The  transition  line  being  a  row  of  grid  points,  the  differ¬ 
ence  equations  for  these  points  degenerate  to  <j>  =  0,  ex- 

AA 

actly  as  does  the  differential  equation.  This,  together 

with  the  side-wail  periodicity  requirements,  yields  <{i  *  const 

on  y  =  0.  From  what  was  said  about  the  n  =  0  case  in  the 

previous  section,  we  are  permitted  to  prescribe  this  constant 

* 

to,  say,  <J>(x,0)  =  0. 

In  this  simple  case,  the  elliptic  and  hyperbolic  regions 

* 

completely  uncouple.  The  difference  equations  in  the  ellip¬ 
tic  region  possess  a  unique  solution,  as  we  shall  show  below. 
The  situation  is  not  much  different  from  the  Laplace  equation, 

*  The  model  is  too  restrictive  compared  to  the  real  nozzle, 
in  this  sense.  The  transition  line  in  the  real  nozzle  is 
almost  a  characteristic,  whereas  in  the  model  it  is  exactly 
a  characteristic. 
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since  the  degeneracy  on  the  x-axis  has  been  properly  taken 
care  of.  Once  the  unique  solution  is  obtained  in  the  ellip¬ 
tic  region,  it  is  trivial  to  show  that  the  equations  in  the 
hyperbolic  region  have  unique  solutions  line  by  line  (stan¬ 
dard  implicit  scheme  for  wave  equation). 

To  prove  unique  existence  for  the  elliptic  region,  let 
the  elliptic  region  have  m  rows  of  interior  grid  points 
in  the  x-direction,  and  n  rows  in  the  y-direction,  fig.  3. 
The  difference  equations  for  this  region  become 


A<j>  =  f 


(14) 


where  the  column  vector  £  denotes  (4*^  <J>12  •  •  •  ^in ’  ^21  ‘  '  ^mn ^ 

and  the  column  vector  f  denotes  the  boundary  data  as  usual. 


The  matrix 


has  the  form 


Q  -I  0  0 

-I  Q  -I  0 
0  -I  Q  -I 


0  -I 


0  0 


(15) 


I- 1  0  0 


..  -I  Q 


where  submatrices  Q  and  I  are  n  x  n  matrices,  and  there 
are  m  x  m  of  them  to  make  up  A. 

For  Neumann  data  on  y  =  -a,  the  matrix  Q  is 


nr+2 


-(n-l)r  2(n-l)r+2  -(n-l)r  0 


-(n-2)r  2(n-2)r+2  -(n-2)r  ..  0 


0  ....  -2r  4r+2  -2r 


(16) 


0  ....  0 


-r  2r+2 


where  r  =  Ax  /  Ay.  For  Dirichlet  data  on  y  =  -a,  Q  is 
the  same  matrix  as  (16),  except  that  the  first  element 
is  replaced  by  2nr+2. 

We  now  use  a  theorem  of  Taussky-Todd  |5j  ,  which  states 

that  for  irreducible  matrices  M,  with  positive  diagonal  ele¬ 
ments.  >0,  and  nonpositive  off-diagonal  elements,  <  0 
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0  for 


for  \f\i  ,  M  is  singular  if  and  only  if  T  M,,, 

]1  A  ** 

every  A  .  Our  matrix  A  for  both  Neumann  and  Dirichlet 
data  on  y  =  -a  satisfies  all  these  conditions,  except 

A^  f  0  for  A  =  n,  2n ,  ...  mn.  Hence  A  is  nonsingu¬ 
lar,  and  the  unique  solvability  of  the  difference  equations 
in  the  elliptic  region  is  proved. 

The  condition  that  y  =  0  be  a  row  of  grid  points  will 
now  be  removed.  To  this  end,  let  the  transition  line  y  =  0 
be  at  a  distance  aAy  from  the  last  elliptic  row  (j  =  n) 
and  a  distance  of  (l-m)Ay  from  the  first  hyperbolic  row 
( j  =  n+1) ,  where  0  <  ot  <•  1  .  Defining  (l-a)4>i  n  +  n+1  * 

substituting  j  =  n  and  j  =  n+1  into  the  first  and  second 
forms  of  (13)  respectively,  and  adding  them  properly,  we  get 

♦t+i  *  *i-i  -  2h  ’  0  • 


Obviously,  is  the  interpolated  value  for  4> ( x , 0 ) ,  and 

this  equation  approximates  4>  =0  on  y  =  0.  This  equation 

should  be  solved  first,  again  with  periodicity  conditions 
on  the  side  walls,  to  give  $ ^  =  const  =  0  as  in  the  pre¬ 
vious  case. 

Then  we  replace  <j>i  n+1  by  -  n>  and  again  have 

just  a  mn  x  mn  system  of  linear  equations  in  the  elliptic 
region.  Matrix  A  once  more  has  the  form  given  in  (15), 
except  that  Q  (for  Neumann  data  on  y  =  -a)  is  now 

(n-l+a)r+2  -(n-l+a)r  0  0....  0 

-(n-2+a)r  2(n-2+a)r+2  -(n-2+a)r  0....  0 


0  0  .  -(l+a)r  2(l+a)r+2  -(l+a)r 

0  0  .  0  -ar  (l+a)r+2 


Again,  the  theorem  of  Taussky-Todd  insures  the  nonsingulari¬ 
ty  of  A.  A  similar  argument  also  holds  for  Dirichlet  data 
on  y  =  -a. 


We  conclude  this  section  with  some  computed  results. 

Fig.  4  shows  a  model  nozzle  for  the  inverted  Tricomi  equa¬ 
tion,  calculated  with  <p  =  0  on  the  transition  line,  periodic 
boundary  conditions  on  the  side  walls,  and  two  cases  of  Neu- 


mann  data  on  y«-a:(i)$»l+  cos  x,  and  (ii)  $y  * 

1  +  cos  x  +  sin  x.  Ax  and  Ay  were  taken  to  be  ti/8  and 
1/5  respectively,  giving  n  *  4  and  m  =  8  in  the  elliptic 
region.  Presented  in  the  figure  are  values  of  <j>  (modeling 
the  axial  velocity)  as  functions  of  x  at  various  values  of 
y.  The  solid  curves  represent  the  analytic  solution  described 

in  Sec.  3,  while  the  circles  are  the  numerical  solution.  The 
accuracy  is  obviously  very  satisfactory. 


Fig.  4.  Calculated  model  nozzle  flow 
compared  with  exact  solution. 


a.  finite  difference  solution  for  the  periodic  tricomi  prob¬ 
lem  (FLOPPY  DISC  MODEL) 


For  the  periodic  problem  of  the  Tricomi  equation,  as 
described  in  Section  2,  we  recall  that  boundary  conditions 
are  needed  at  both  y  *  -a  and  y  *  +b,  and  that  a  unique 
solution  exists  only  when  the  domain  does  not  admit  eigen- 
solutions  corresponding  to  wave  numbers  n  ^  N.  An  obvious 
method  to  solve  such  problems  is  indeed  to  solve  equation  (6) 
as  an  ordinary  differential  equation,  and  treat  the  eigen¬ 
functions  as  part  of  that  problem.  Here  we  propose,  however, 
an  alternate  procedure,  starting  directly  with  a  finite 
difference  grid  for  the  partial  differential  equation. 

We  denote  by  $  (x)  and  4>0(x)  the  boundary  values  given 
on  y  =  -a  and  y  =  +b  respectively,  and  we  denote  by  ^(x) 
the  function  $(x,0)  on  the  transition  line.  If  were 

known,  the  solution  of  the  difference  equations  in  the  ellip¬ 
tic  region  will  exist  uniquely,  since  the  matrix  A  will  be 
nonsingular  for  exactly  the  same  reasons  as  with  the  inverted 
Tricomi  equation  discussed  in  the  previous  section.  Then  we 
can  apply  the  Murman-Cole  scheme  in  the  hyperbolic  region  as 
before,  and  obtain  the  entire  solution  up  to  <j)(x,b),  which 
should  be  equal  to  ^(x).  But  4>(x,b)  is  a  linear  trans¬ 
formation  of  4>0(*)  and  of  <J>^(x)  summed,  which  gives  pre¬ 
cisely  a  condition  to  determine  01(x): 

<(>(x,b)  =  S(t>0  +  T^  =  4>  2 

Or, 

*1  =  T  1(  ^2  "  SV  (17) 

Here  S  and  T  are  the  matrices  representing  the  linear 
transformations,  and  the  method  works  only  if  T  has  an  in¬ 
verse.  This  condition  is  exactly  the  finite  difference  coun¬ 
terpart  of  the  nonoccurrence  of  the  eigenfunctions. 

In  actual  computations,  we  first  solve  the  entire  prob¬ 
lem  with  ^(x)  =  0,  and  get  4>(x,b)  =  S<f>Q.  We  then  solve  a 
series  of  problems,  using  4>0  =  0 .  and 

4,i(i)  =  (0,0,  .  .  ,1,0,  .  .0)  i  =  1  , - n 

(i.e.,  $  has  1  ln  tho  i-th  Placc*  and  0  elsewhere). 
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T  is  constructed  in  a  simple  manner  from  the  solutions  of 
these  problems  at  y  =  b.  The  nonsingularity  of  T  is  pre¬ 
cisely  the  criterion  for  the  unique  solvability  of  the  entire 
problem. 

Using  this  procedure,  we  calculated  an  example  using  <t>(x, 
b)«  <j>2 (x )  *  1  +  sin  x  and  <j>(x,-a)  *  <t>Q(x )  =  1+  cos  x.  The 
results  are  shown  in  Fig.  5.  Again,  the  solid  curves  are  the 
exact  solutions  from  Sec.  2,  and  the  circles  are  the  computed 
results.  The  accuracy  is  again  satisfactory. 

Whether  or  not  this  procedure  is  more  efficient  than  solv¬ 
ing  the  ordinary  differential  equations  is  debatable,  but  this 
procedure  is  interesting  in  its  own  right,  and  permits  a  com¬ 
pletely  different  way  to  look  at  the  problem. 


-'it  -n/2  0  ir/2  it  x 
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